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We explore the quantum dynamics of a one-dimensional trapped ultracold ensemble of bosonic 
atoms triggered by the sudden creation of a single ion. The numerical simulations are performed by 
means of the ab initio multiconfiguration time-dependent Hartree method for bosons which takes 
into account all correlations. The dynamics is analyzed via a cluster expansion approach, adapted 
to bosonic systems of fixed particle number, which provides a comprehensive understanding of the 
occurring many-body processes. After a transient during which the atomic ensemble separates 
into fractions which are unbound and bound with respect to the ion, we observe an oscillation 
in the atomic density which we attribute to the additional length and energy scale induced by the 
attractive long-range atom-ion interaction. This oscillation is shown to be the main source of spatial 
coherence and population transfer between the bound and the unbound atomic fraction. Moreover, 
the dynamics exhibits collapse and revival behavior caused by the dynamical build-up of two-particle 
correlations demonstrating that a beyond mean-field description is indispensable. 
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I. INTRODUCTION 

The theoretical description of degenerate atomic quan¬ 
tum gases has advanced significantly in the past two 
decades. It relies particularly on the separation of 
length scales of the external trapping potential, the inter¬ 
particle distances and the range of atomic interactions. 
In the ultracold regime, this allows to approximate the in¬ 
teraction by a contact potential (or simply pseudopoten¬ 
tial) [l], which has proven tremendously powerful when 
applied to bosonic gases. The latter holds especially for 
the mean-field approximation leadingto the well-known 
Gross-Pitaevsikii (GP) equation 0 [3] which describes 
the condensed state of a weakly interacting bosonic gas 
yielding a variety of phenomena such as collective exci¬ 
tations, solitons, and vortices 3,0. For state-of-the- 
art methods, such as the density-matrix renormaliza¬ 
tion group 0 and the multiconfiguration time-dependent 
Hartree method for bosons (MCTDHB) Jm as well as its 
multi-layer extension (ML-MCTDHB) 00, the pseu¬ 
dopotential represents an essential simplification for the 
quantum dynamical description of many-boson systems 
in order to investigate the physics beyond the mean-field 
approximation. Even though the short-range interaction 
has proven to lead to exotic states of quantum matter like 
supersolidity [Tj ) or crystalline phases [lJL] and therefore 
represents an important case, not all ultracold atomic 
systems exhibit these short-range interactions. An ex¬ 
ample for the latter are chromium atoms which possess 
long-range magnetic dipolar interactions mm- 

Hybrid atom-ion systems represent a specific class of 
systems with a new scale of interactions due to the inter- 
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play between the charge of an ion and the induced dipole 
moment of a neutral atom [l4j]. They have recently be¬ 
come available experimentally [ 1 ri 1 1 81 and attracted in¬ 
creasing interest. Most of the current experiments, how¬ 
ever, are based on the Paul trap scheme, whose draw¬ 
back is the so-called micromotion which so far prevents 
from reaching the ultracold regime as shown in theoreti¬ 
cal classical and quantum analyses pjl [2C|. Correspond¬ 
ing studies showed that a large ion-atom mass ratio might 
help circumventing this limitation fl9l42ll | , although it is 
not yet clear whether the s-wave regime can be reached. 
Still, a recent detailed study has shown that specifically 
for the atom-ion pair Li-Yb + the ultracold regime can 
be reached experimentally (22]. Given these findings, it 
is desirable to extend the theoretical understanding of 
ultracold neutral quantum many-body systems by the 
presence of ions. This attractive interaction induces, ad¬ 
ditional to the confinement, a further length and energy 
scale and it is therefore expected to lead to intriguing ef¬ 
fects such as the formation of molecular ions 23| and ion- 
induced density bubbles in the atomic cloud [241 ] . Apart 
from this fundamental point of view, these hybrid sys¬ 
tems show versatile applications in quantum information 
processing as, for example, the controlled creation of en¬ 
tanglement [2l1. [25l] . the realization of quantum gates [Hj, 
and the simulations of solid-state systems f27| . 

In a previous study, we have investigated in detail the 
ground-state properties of an atom-ion hybrid system, 
consisting of a single static ion in the center of a bosonic 
atomic cloud. The dependence of relevant observables 
on the atom number and the interaction strength has 
been analyzed and we showed that the presence of an 
ion strongly affects them [28|. For weakly interacting 
atoms, we found that the ion impedes the transition to 
the Thomas-Fermi regime while for strong atom-atom in¬ 
teraction it modifies the fragmentation behavior depend- 
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ing on the atom number parity. In the present work, we 
explore the impact of the second length and energy scale 
generated by the atom-ion interaction onto the dynam¬ 
ics of the atomic cloud. In particular, we envision the 
scenario in which we create a single ionic impurity in an 
ensemble of N interacting atomic bosons. If the atomic 
cloud, initially prepared in the ground state of a harmonic 
trap, and the ion are suddenly brought into contact, this 
process resembles, to some extent, the sudden ionization 
of a single impurity atom within the atomic cloud. Here, 
however, we still neglect the ionic motion which is jus¬ 
tified in case the ion is tightly trapped. We will see in 
the following that the second length and energy scale in¬ 
duces a coherent oscillation between states bound in the 
atom-ion potential and states of the harmonic trap. This 
oscillation occurs in addition to the usual harmonic exci¬ 
tation and reveals a collapse and revival behavior caused 
by the dynamical build-up of correlations. 

This work is organized as follows. In Sec. IE we define 
the setup and explain the model of our hybrid atom-ion 
system. In addition, we introduce a cluster-expansion 
scheme which enables us to distinguish the single-particle 
dominated physics from the processes induced by the 
build-up of (quantum) correlations. In Sec. Mil we ex¬ 
plore the dynamical evolution of the hybrid system and 
identify the most important excited modes by means of 
our cluster-expansion approach. Sec. he contains a brief 
convergence analysis which shows the necessity to use an 
advanced method like MCTDHB. Finally, we summarize 
our findings in Sec. IE together with the conclusions and 
an outlook on future investigations. 


II. MODEL AND THEORETICAL APPROACH 

In the following, we introduce our model and define 
the process to initiate the dynamics. Further, we briefly 
outline the MCTDHB used for the simulations, define 
important quantities for the analysis, and introduce our 
cluster expansion approach which we exploit in order to 
analyze the many-body wave function. 


A. Model of the System 

We consider N interacting bosonic atoms in a one¬ 
dimensional (ID) harmonic trap at zero temperature ini¬ 
tially prepared in the ground state of the system which we 
compute by imaginary time propagation |29| of an initial 
guess wave function. The short-range intra-atomic in¬ 
teraction is modeled by a contact pseudopotential. Into 
this atomic cloud, we immerse a single trapped impurity 
atom which does not interact with the other atoms. This 
could be achieved by tuning the inter-atomic interaction 
to zero by exploiting a Feshbach resonance [3(j • At time 
to = 0 this impurity atom is ionized by a laser pulse such 
that a single ion is created within the atomic cloud. The 
ionizing laser pulse is assumed to be far detuned from 


any possible resonance of the atomic cloud. Moreover, we 
assume pulse durations such that the ionization process 
takes place on much faster time scales than any possi¬ 
ble response of the atomic cloud. This allows us to treat 
the ionization as an effectively instantaneous process. If 
the atomic impurity is trapped in a sufficiently deep and 
tight trap, the proposed scenario can be indeed achieved 
experimentally, as recently reported for optical trapping 
of a single ion [3l|]. Thus, at time to = 0, the ion is as¬ 
sumed to be statically trapped at z\ = 0 in the atomic 
cloud. The motional excitations following this ionization 
process result in a rich dynamics which we analyze in Sec. 

HE 

The interaction between the ion at z\ and an atom at 
z\ behaves in ID at large distances as —ae 2 /(2 (za — Zi) 4 ) 
up to a minimal cutoff distance Rid f32j, where a is the 
polarizability of the atoms and e the elementary charge. 
For numerical many-body simulations with MCTDHB it 
is more convenient, however, to define a model potential 
as |28| : 

W^voe- 2 --^. (1) 

Here z denotes the relative coordinate z\ — z\. The model 
parameters vo, 7 , and uj are determined by the short- 
range quantum defect parameters. The above potential 
asymptotically approach the — 1/z 4 behavior at large dis¬ 
tances, whereas at short distances it has a barrier which 
is designed such that the quantum defect theory results 
for the atom-ion scattering [33| can be reproduced. Ad¬ 
ditionally to the harmonic confinement, this interaction 
introduces a second length R* = \Jae 2 m/h 2 and energy 
scale E* = h 2 /(2mR* 2 ) to the system, where m is the 
mass of a single neutral atom. 

In summary, the system can be described by the fol¬ 
lowing many-body Hamiltonian (in E* and R* units) 
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with the harmonic trap of frequency wo and characteristic 
length l = yJhJijrujJo) / R *, and the intra-atomic contact 
interaction strength g. The ionic part is switched on at 
time to by a step-function 6{t — to)- Thereby, we term the 
Hamiltonian of a single boson as R. In the following, we 
denote the stationary eigenenergies of hi(t > to) for fixed 
t : with ej and the corresponding single-particle eigen¬ 
states with (/>j(zi) such that fn(t > to)<jPj(zi) = 
Hereafter, for the sake of numerical convenience, we set 
l = 0.5 R*, which would corresponds to uiq = 2n ■ 3.3 kHz 
and l = 188 nm for s 'Rb atoms. We will see that this 
choice does not affect qualitatively the observed dynam¬ 
ics. Furthermore, we choose a weak intra-atomic inter¬ 
action strength g = 2 E*R* and small atomic ensembles 
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FIG. 1. The effective potential (gray shaded area in back¬ 
ground) consisting of the harmonic trap and the atom-ion in¬ 
teraction potential for the 87 Rb atom. Further, we show the 
single-particle energy levels tj indicated by straight horizon¬ 
tal lines. In addition, the energetically lowest single-particle 
state 0? (z) (dark gray area) and the trap state 0s(a) (light 
gray area) are sketched with arbitrary but equal scaling. Pos¬ 
sible types of dynamics occurring in the system are indicated 
by arrows. 


consisting of N = 2 up to N = 10 neutral atoms. Be¬ 
sides, we fix the model parameters to to = 80(i?*) -4 , 
v 0 = 3 u) and 7 = 4\/10 w. We refer here to Ref. 0 for 
a detailed discussion of the chosen parameters as well as 
for the experimental conditions needed for the quasi-ID 
regime. Our above choice leads to two bound states for 
the atoms in the atom-ion potential which are localized 
on both sides of the ion but vanish at Z\. Even though 
we neglect the motion of the ion, we term the two states 
below E = 0 bound states while we refer to the remain¬ 
ing states as trap states (with E > Hujq/2). In Fig. [I] 
we show the total potential for the atoms together with 
the energies 6j as well as the lowest single-particle state 
bound in the atom-ion potential, and the state 
0° ( 2 ). The arrows indicate the possible processes that 
can occur: within the trap (green arrow), in the atom- 
ion potential (white arrow), and between those two scales 
(oblique magenta arrows). 


B. Theoretical Approach 


We explore the quantum dynamics of the many-body 
system described by the wave function | SR) by means of 
the numerically exact ab initio method MCTDHB. Its 
main idea is that by using m time-dependent variation- 
ally optimized single-particle basis functions, the number 
of basis functions can be kept rather small. More pre¬ 
cisely, the many-body wave function |4/) for N bosons is 


expanded in bosonic number states |n(t)) 

|*(i)) = 5> n (*)|n(i)> (3) 

n\N 


in order to take into account the indistinguishability 
of the bosons. Note that in a number state |n(f)}, 
each boson occupies one of the m time-dependent single¬ 
particle functions (SPFs) |'F 3 (i)} and that the vector 
n = (ni,--- ,n m ) contains the occupation numbers nj 
of every SPF. Besides, the sum in Eq. d3]) goes over all 
possible n with JT n 7 = N , which is denoted by the 
symbol n|IV. With this ansatz for the many-body wave 
function, the temporal evolution of the wave function 
|4/(t)) is obtained by means of the Dirac-Frenkel varia¬ 
tional principle [34L |35| which guarantees a variational 
optimal many-body solution. We would like to empha¬ 
size that not only the coefficients A n (t) but also the SPFs 
|T 7 (f)) are adapted in time to the many-body dynamics 
in order to allow for the largest possible overlap between 
the ansatz Q and the true many-body wave function. 
We refer for a detailed description of the method to Refs. 

0-S- 

The analysis of the full many-body wave function |^) 
is generally a complicated task due to the underlying 
high dimensionality resulting from the N degrees of free¬ 
dom. In order to analyze the time-dependent many-body 
wave function in detail, we inspect two different quanti¬ 
ties: The one- and two-particle reduced density matrices 
which allow for the investigation of spatial coherence and 
correlations [36|, and clusters enabling us to analyze co¬ 
herence and correlations in terms of any single-particle 
basis. We transfer and adapt the notion of clusters from 
Ref. [37j to bosonic systems of few particles. 


1. Density Matrices 

The reduced one- and two-particle density matrices are 
defined via the expectation values of the field operators 
'l't(x) and d* (x) as p\(x,y,t) = (4d(x,f)'k {y,t)) and 
p 2 (x,y,y',x',t) = {&(x,t)&(y,t)$ (y',t)T (. x',t )), re¬ 
spectively. Their spectral decomposition can be written 
in terms of the natural populations A j (t) and the natural 
orbitals t) 

P(x,x',t) = ^Aj(t)$j(x,t)$j(x , ,t), (4) 

i 

and the natural populations 7 j(t) and natural geminals 
<Pj(x, x', t) 


P 2 (x,y,y',x',t) = 5^7 j(t)$*(x,y,t)$j(x',y',t), (5) 

3 

respectively. This natural representation of the reduced 
density matrices has several advantages. For example, 
the natural populations can be used to identify the de¬ 
gree of fragmentation of the system [38J] and to judge the 
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convergence of our numerical simulations [39j . Further, 
the natural orbitals build always a very suitable single¬ 
particle basis set for the description of the dynamical sys¬ 
tem, hence in many cases only a few basis functions are 
needed to represent the (many-body) wave function in 
this basis. Despite these advantages, the analysis and a 
deep physical understanding of the quantum dynamics in 
this basis is very difficult, since the natural orbitals are 
time-dependent eigenfunctions of the density matrix in 
its spatial representation. Given this, we use instead the 
so-called clusters for the detailed analysis of the many- 
body dynamics. We define in the following the notion 
of clusters in the framework of the cluster-expansion ap¬ 
proach. 


2. Definition of Clusters and Correlations 

The cluster-expansion approach is a very powerful 
technique to describe the quantum dynamics of in¬ 
teracting many-body systems because it allows for a 
systematic truncation of the Bogolyubov-Born-Green- 
Kirkwood-Yvon (BBGKY) hierarchy [4(j. By expanding 
M-particle expectation values into cumulants or corre¬ 
lated clusters, a consistent theory up to the sing le-, two-, 
or even M-particle level can be developed [4J|. In this 
spirit, we shall use the cluster-expansion approach in or¬ 
der to separate the many-body dynamics, obtained by 
means of MCTDHB, into single- and two-particle contri¬ 
butions which will enable us to get more physical insight. 
In particular, it will be useful for the identification and 
classification of the most relevant excitations of the sys¬ 
tem that are created during the dynamical evolution (see 
Sec. m 

To this end, we will first briefly review the cluster- 
expansion approach, which will be helpful for a better 
understanding of our modifications to the traditional ap¬ 
proach. Indeed, we shall adapt the traditional cluster- 
expansion to bosonic systems with a fixed particle num¬ 
ber. This constraint is particularly important for atomic 
ensembles of only a few particles. 

To begin with, let us consider an arbitrary orthonor¬ 
mal time-independent single-particle basis <pj (x) and the 
annihilation and creation operators aj and a ] for the cor¬ 
responding single-particle state, respectively. Then, the 
expectation value of any observable can be expressed in 
terms of the M-particle expectation values or M-particle 
clusters 

M M 

<*|Af Jik |*> = <Mj, k > = 

i i 

with the index sets j = { ji , ...,j’m} and k = {/q,..., Izm} 
and M < N. These are nothing else but the matrix el¬ 
ements of the M-particle reduced density matrix in an 
arbitrary basis and can be, most conveniently for MCT¬ 
DHB, obtained from the spectral representation of the 


M-particle reduced density matrix, as shown for M = 1 
and M = 2 in App. [A]. 

Now, the single-particle properties of the system are 
given by the single-particle clusters (at a-), also called 
singlets. We distinguish between occupations fj = 
(aja ■), which describe the population of the state 

and coherences pij = (i ^ j) which can be under¬ 

stood as a transition amplitude between the i- th and the 
j-th state. 

Starting from the singlets, the cluster-expansion is re¬ 
cursively build-up based on the consistent factorization 
of an M-particle cluster into independent particles (sin¬ 
glets), correlated pairs, correlated three particle clusters, 
up to correlated M-particle clusters [37l. [42^. that is: 

(1) = [<l)]s (7) 

(2) = [(2)]s + A(2) c (8) 

(3) = [{3)]s + [A(2) c (1)]d + A(3)c (9) 


Here the terms [(M)]s represent the single-particle con¬ 
tributions, while the terms A(M)c contain the corre¬ 
lated part of the M-particle cluster. Note that we omit¬ 
ted the indices for brevity such that all cluster prod¬ 
ucts, denoted in square brackets, include a sum over all 
unique permutations. For instance, for the two-particle 
clusters (a|,aja ? ,a fe ,), the single-particle contributions are 
defined as [<2)] s = [{ala\a q ,a k ,)} s := {a\a q ,){a\a k ,) + 
(a£a fc/ )(aja /). Given this, the correlated part of the two- 
particle cluster is given by A(2)c := (2) — [(2)]s- Now 
the correlated clusters A(M)c with M > 2 can be de¬ 
termined recursively which completes the formulation of 
the cluster expansion. 

At this point a remark is in order: For bosonic 
systems with a fixed number of particles (i.e., in a 
number-conserving theory), even in a GP type mean- 
field state (i.e., with only one single-particle orbital), the 
term A(M)c is non-zero because of the bosonic symme¬ 
try. This implies that the systematic truncation of the 
BBGKY hierarchy, for which the correlated parts of any 
M-particle cluster beyond a certain size (M > Mt) have 
to be neglected, cannot be performed in this case. Note 
that this issue does not arise neither for fermions nor for 
bosons with particle number fluctuations. Thus, in order 
to circumvent this problem, we shall introduce a slightly 
different definition of the correlated parts A(M)c- Our 
strategy will be to define the correlated parts of any M- 
particle cluster in such a way that all the (M— l)-particle 
contributions are indeed removed. As a result of such a 
strategy, for example, the correlated parts automatically 
vanish in any mean-field state. 

To this end, let us note that an M-particle cluster con¬ 
tains all the information about the M'-particle clusters 
with M' < M. This can be easily seen by using the 
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recursive relation between the density matrices 

PM—! ( x 1j ' t X M—li x M — 1’*'* ?*^l) = 

d X M PM (. X 1^ ; X M ; X M ; M — 15 " ? *^i) 

( 10 ) 

for M > 2, which leads to 

« M “ ‘W = JV-M+1 (U) 

9 

In order to identify the correlated part of an M-particle 
cluster, we decompose the M-particle cluster into two 
parts: one consisting of all contributions from clusters 
with M' < M , denoted by (Myk)<M, and one which 
contains only the M-particle contributions, that is, 

(Mj, k ) := (Mj. k )<M + A(M jik ). (12) 

There are several ways one could perform such a 
decomposition. For instance, in the traditional 

cluster-expansion approach, as we have discussed 
above, one would choose the following definition for 
the single-particle part of the two-particle cluster: 
( 2 {k, q }{q',k '})<2 = (d{a q ,){d\a k ,) + {a\a k ,){a\d q ,). Here, 
however, we shall define the term (Mj ik )<M by requiring 
that it has to fulfill the condition 

<(M - l)j,k) = N _ l M + 1 5Z^ M fj.9},{9,k})<M- (13) 

9 

This expression is assumed to hold for any many-body 
quantum state and implies that = 0 

pH! . Besides, if we consider, for instance, the case M = 2, 
then we see that the right-hand side of Eq. m accounts 
only for single-particle contributions of two-particle clus¬ 
ters. 

In order to find a definition for (My k )<M such that Eq. 
(fl3l) is fulfilled, let us first investigate, as an example, the 
case where |\&) is a general mean-field state, which is de¬ 
fined as a single permanent. More precisely, given some 
single-particle orbitals x.j ( x ) with the associated creation 
operators dj, a mean-held state is defined as a single per¬ 
manent like |MF k ) = nP Ckj) l vac ) with |vac) being the 
vacuum and k = {Aq, ...,fcj\r} (for the commonly known 
GP state ki = AqVi). For such a state the two-particle 
clusters can be written as 

(ala\a q ,a k ,) M F =(d\a q ,)MF(d\d k ,)M¥ 

+ (PU fc /)MF (a\a ? ,)mf 
+ A^ f (k, q, q , k 1 ) (14) 

with the bosonic correlations [44t | given by 
A b F (k, q , q', k') = - ^{xMk){Xj\4>q) 

3 

{(l>k'\Xj}(<t>q'\Xj)nj{nj + 1). (15) 


Here nj is the number of bosons in the single-particle 
state Xj( x ) defined via rij = ’Yh i Sk i ,j {$i,j is the 
Kronecker-Delta). It is easy to verify that for such 
a mean-held state Eq. m is fulfilled if we set 
{a\a\a q ,a k ,) < M = (aj[dta g ,a fc ,)MF- With this choice, the 
term A(2^ kq y,{q',k'}) is zero in a mean-held state per 
dehnition. In this spirit, the terms A(Mj. k ) can be un¬ 
derstood as the M-particle correlations. 

Now for a general many-body quantum state we re¬ 
place in Eq. (11511 the mean-held occupations nj and the 
mean-held basis functions Xj( x ) by the natural popu¬ 
lations A j and the natural orbitals &j{x), respectively, 
yielding the following analogue expression: 

A b {k,q,q’,k') = - 

3 

<^;>< 09 «A,(A,+i). (16) 

Note that in general A j is not an integer. Hence, we 
dehne the single-particle contributions of a two-particle 
cluster as 

(aW q a q ,a k ,) <M :={a\a q ,){a\a k ,) + (ala k ,){a\a q ,) 

+ A B (fc,<?,g',fc')> (17) 

whereas the two-particle correlations, also called dou¬ 
blets , are dehned as 

Qkqq’k' =A(a\a\a q ,a k ,) 

■.={a\d\a q ,d k ,) 

- {a\a q ,){a\a k ,) - {a\a k ,){a\a q ,) 

-A B (k,q,q',k'). (18) 

Here some considerations are in order. Our choice for 
the non-correlated part of the two-particle cluster [Eq. 
m] indeed fulfills the condition (ITTIh which justihes the 
above replacement of the mean-held occupations and ba¬ 
sis functions by the natural populations and natural or¬ 
bitals, respectively. Although we focused here on the two- 
particle clusters, we note that expressions like Eq. m 
can be obtained for any M-particle clusters, too. How¬ 
ever, for the present study the singlets and doublets are 
sufficient to understand the dynamics of the system. Fur¬ 
ther, we would like to stress that the decomposition into 
singlets and M-particle correlations does not correspond 
to a separation into mean-held and beyond mean-held 
contributions. But even if the condition m only sep¬ 
arates the clusters into singlets, doublets, three-particle 
correlations, etc., it enables us to identify genuine correla¬ 
tions of any many-body quantum state which are beyond 
mean-held. Indeed, in the limit of a mean-held state Eq. 
m boils down to Eq. m , and therefore all correlated 
parts A(2 {k, q }.{ q ’,k'}) = 0 vanish. This would not be 
possible with the traditional cluster-expansion approach. 

Finally, we would like to highlight that one could in¬ 
stead search for the best mean-held state [45] and then 
separate the M-particle clusters into mean-held part and 
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contributions beyond that. It turns out, however, that 
such a choice does not satisfy Eq. m- This shows that a 
mean-held approximation of a real many-body state does 
not necessarily result in the exact single-particle proper¬ 
ties of the system. 

3. Singlet Dynamics 

Even though we are able to derive the time evolution 
of every M-particle cluster from our MCTDHB solutions, 
it is worth to investigate the equations of motion of the 
clusters since they make it possible to understand the 
coupling between the singlets themselves and between 
singlets and doublets. The dynamics of the M-particle 
clusters can be derived from the Heisenberg equation of 
motion and the above outlined definitions. For the anal¬ 
ysis of the dynamics in Sec. m we focus onto the equa¬ 
tions of motion of the singlets. One can easily show that 
the time evolution of the coherences is given by 

= (O “ *i)Pii + s jiifi - fj ) 

+ E \pjqPiq ~ ^ iqPqj\ 

q=£i,j 

+ T® +Tjj (19) 

while the one of the populations is governed by 

=2 ^2 Im [T, iq p iq ] + Im [T? + T^] . (20) 

q^i 

Here we used the definition of Eq. m for the two- 
particle cluster that appears in the corresponding Heisen¬ 
berg equation, because of the atom-atom interaction and 
the asterisk for the complex conjugation. Besides, in the 
above equations, we have introduced the renormalized 
single-particle energies e, = ei + E^, the singlet couplings 
E ij = 2 Ylkq Vfeyq(aJ.d g ), the coupling to the bosonic cor¬ 
relations 

r® = Wkjqk' A B (b k, q, k') - Vkqik' Ab(<?, k,j , fc')], 

kk' q 

( 21 ) 

and to the doublets 

r ij — E \Ykjqk' Qikqk' Vkqik' Qqkjk'\ •> (22) 

kk' q 

and the interaction matrix elements 

Vijj'i' = 9 j 4>i{x)<l>*j{x)^j>{x)<l>i>{x)Ax. (23) 

Note that the equations of motion for the singlets form 
a system of non-linear differential equations, since the 
mean-field couplings are also defined by the singlets 
themselves. Furthermore, they are coupled to the dou¬ 
blets (via T^) which can be interpreted as a source term 


(or inhomogeneity) for the singlets. On the other hand, 
the dynamics of the doublets depends on the singlets as 
well as on the three-particle correlations A(3), which ren¬ 
ders the coupling time-dependent. We would like to em¬ 
phasize that the coupling to A(3) is a manifestation of 
the BBGKY-hierarchy. 

For the truncation of the hierarchy at the singlet level, 
one has to ensure that the contribution of the doublets 
to the singlet dynamics remains very small during the 
time interval of interest, which would imply that one can 
set T^ = 0. In addition, we remind here again that 
in order to obtain a closed singlet theory, the coupling 
to the bosonic correlations has to be included which is 
the price to pay for our choice of factorization. In the 
following, we show how a consistent and closed singlet 
theory can be accomplished. At first, by means of Eqs. 
(1131) and (ED- we obtain an exact relation among the 
bosonic correlations and the singlets of the system: 

E AB ( fc ’ Q’ “ E + (Pla q )(a\a k ,) . 

Q Q 

(24) 

Then, by using this relation and by noticing that 
Ab (i,k,q,k r ) ^ 0 only if all indices are even or if all 
are odd [see Eq. m and note that the natural orbitals 
have a defined parity in our setup] we can approximate 


as 



bAE 

Q 

Vqqqi 

\ k 

Vjqqq ^ 



With this expression for the coupling to the bosonic cor¬ 
relations, the equations of motion of the singlets can 
be completely decoupled from higher than single-particle 
clusters such that a consistent and closed singlet theory 
is obtained. In the following, we will see the power of 
such a singlet theory in the analysis of the complicated 
many-body dynamics, especially in combination with the 
MCTDHB method. 

III. DYNAMICAL EVOLUTION 

We investigate now the dynamics induced by the in¬ 
stantaneous creation of an ion in an atomic cloud. First, 
we analyze the one-body density (matrix) and the com¬ 
ponents of the energy as well as the many-body excita¬ 
tion spectrum. The observations are then discussed and 
analyzed in detail in terms of our singlet-doublet theory 
developed above. 

A. Observations 

In Fig. 0 we show exemplarily the temporal evolu¬ 
tion of the one-body density p(z,t ) = p(z,z,t ) of the 
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t (units of h/E*) 


FIG. 2. Time evolution of the one-particle density p(z, t) and 
the components of the total energy per particle for a system 
with N — 2. The green, cyan, and magenta lines represent the 
trapping, kinetic, and ionic energy per particle, respectively. 


atomic cloud for g = 2 E*R* and N = 2. In order 
to give a feeling of the actual time scales involved, we 
note that for 87 Rb we have h/E* ps 0.39 ms, while for 
7 Li it corresponds to h/E* ss 0.001ms. For short times, 
the suddenly created ion captures very quickly most of 
the atomic cloud in its bound states. The remaining 
atomic density fraction is emitted as a beam into the 
outer region of the harmonic trap. Consequently, this 
fraction is decelerated (t < 0.2 h/E*) and back reflected 
(0.2 h/E* < t < QAh/E*) by the harmonic confine¬ 
ment. Subsequently, this sequence repeats with approx- 
imatively constant frequency. Furthermore, a second 
faster oscillation in the density fraction captured within 
the ionic potential is visible (see the holes in the density 
plot at z « ±R*/ 2). Additionally, we show in Fig. [2] 
the components of the energy per particle. The trapping 
energy (green line) perfectly oscillates with a single fre¬ 
quency which coincides with the oscillation frequency of 
the outer density fraction. In contrast, the ionic energy 
(magenta line), which represents the expectation value 
of the ionic potential ©, oscillates with a higher fre¬ 
quency matching the inner density oscillation. On top of 
this oscillation, we observe a short pulse when the outer 
fraction “crashes” into the inner density part. Note that 
these events are not visible in the trapping energy since 
the harmonic trap energy is negligible in the vicinity of 
the trap center. The kinetic energy (cyan line) can be un¬ 
derstood as the negativ sum of trapping and ionic energy, 
since the interaction energy (not shown) is comparably 
small and the total energy is conserved during the dy¬ 
namics. Hereafter, we term the oscillation of the outer 
fraction of the atomic cloud harmonic and the inner frac¬ 
tion ionic oscillation. Note that the above observations 
are qualitatively independent of the atom number N such 
that Fig. [5] is representative also for larger N. 

After multiple oscillation periods (see Fig. [31), we ob¬ 
serve that the ionic energy (magenta line), and therefore 
the ionic oscillation, exhibits a clear collapse and revival 
behavior on this long time scale. On the other hand, the 
trapping energy (green line), and thus the harmonic os¬ 
cillation, becomes only slightly damped. While for N = 2 
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FIG. 3. Time evolution of the various components of the 
energy per particle for N = 10. The red, green, cyan, and 
magenta lines represent the interaction, trapping, kinetic, and 
ionic energy per particle, respectively. 

these two aspects of the long-time behavior can be barely 
seen, it becomes strongly pronounced for larger particle 
numbers. 

The collapse and revival behavior can also be observed 
in the long time dynamics of the atomic density. Figure 
[4] shows the disappearance of the ionic oscillation during 
the time interval [1.5, 3.5 ]h/E* and its recurrence around 
t ss 4.0 h/E*. This effect seems to be directly connected 
to the loss and regain of spatial coherence between the 
inner and outer density fraction which can be observed in 
the snapshots of the one-particle density matrix at times 
t = 0.54 h/E* (left panel), t = 2.20 h/E* (middle panel), 
and t = 4.61 h/E* (right panel) in Fig. [I] Below, we 
will understand the relation between the ionic oscillation 
and the spatial coherence in detail through the singlet- 
doublet analysis. 

Let us now discuss the frequencies that are involved 
in the dynamics. To this aim, we investigate the fidelity 
defined by the overlap of the many-body wavefunction at 
time to = 0 with the one at time t j46j : 

= ( 26 ) 

Since this fidelity F(t) can be understood as the expec¬ 
tation value of the time-evolution operator, thus an N- 
body operator, its Fourier transform contains informa¬ 
tion of all involved excited eigenstates of the interact¬ 
ing N atom system. Due to the discrete nature of the 
spectrum and because in our numerical simulations the 
propagation time T is finite, it turns out to be more ef¬ 
ficient for the computation of the Fourier transform to 
use the compressed sensing (CS) method H3, With 
compressed sensing, one can indeed obtain a resolution 
in frequency space better than Aw = 4S- [49j. To this 
end, we have used the matlab package “SPGL1” for com¬ 
pressed sensing from Ref. [50(. The alg orithms used in 
this package can be found in Refs. [5l], [52|] . In Fig. El 
we show the Fourier transform F(ui) of the fidelity (red 
continuous line). Several prominent resonances become 
apparent. We observe one dominant mode at a frequency 
u> ss 34 E*/h. This corresponds to the ionic oscillation 
frequency. In contrast, the harmonic frequency can not 
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FIG. 4. (Top panel) Time evolution of the one-particle density p(z, t ) for N = 10. Snapshots of the one-particle reduced density 
matrix p{z,z',t) (lower panels) at times t = 0.54 h/E* (left panel), t = 2.20 h/E* (middle panel), and t = 4.61 h/E* (right 
panel) which are indicated in the main top panel as white vertical lines {h/E* = 0.39 ms for 87 Rb atoms). 


directly be found in the spectrum. Furthermore, we see 
that modes with high frequency are excited and seem to 
have an equidistant spacing. In addition to this, there is 
a low energy mode around u « 15 E*/h, whose origin we 
will explain in the subsequent sections. 


B. Singlet Dynamics 

Let us start the analysis of the many-body spectrum 
in terms of the singlets. Therefore, we choose from here 
on the non-interacting single-particle functions as 

the basis (f>j{x ) in which the clusters are expressed. In 
general, the identification of the modes corresponding to 
the observed resonances is a very complicated task. Nev¬ 
ertheless, we were able to identify the most important 
modes by linearizing the equation of motion for the sin¬ 
glets [see Eqs. (flUl) and (O0l) . and App. [B]. The obtained 
energies together with the “exact” spectrum belonging 
to the Fourier transform of the fidelity [see Eq. ([26]) ] 
are shown in Fig. [5] (blue circles). One observes a very 
good agreement of the obtained peak positions with the 
Fourier spectrum of F(t) (red continuous line). On the 
other hand, the relative heights of the peaks can be only 
obtained approximatively by our singlet theory (see App. 
El) showing qualitative agreement with the Fourier spec¬ 
trum. Further, we can identify which singlet has the 
dominant contribution to a certain resonance (see App. 
[B]for details). This dominant singlet is indicated in Fig. 
[5] at the corresponding peak. 

We find that the most dominant frequency corresponds 
to p \3 = (a\a 3 ) excitations. These are coherences be¬ 
tween the number states \N, 0, 0,...) and |iV— 1, 0,1, 0,...) 
(in the basis of the non-interacting single-particle states 
(,b°j{z)) which oscillates with the frequency uq = (£3 — 
6 i )/h = 38.7 E*/h for g = 0. Thus, the ionic oscillation 
corresponds to an oscillation between a bound state and 
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FIG. 5. (Main panel) Excitation spectrum for N = 5 and 
g = 2E*R*. The red line represents the actual Fourier spec¬ 
trum of the fidelity F(t). The blue circles mark the single¬ 
particle energies (see text). We indicate the corresponding 
dominant singlet pk q at each peak by the label pij for the 
sake of better readability. The dashed vertical lines illustrate 
the non-interacting limit, thus correspond to g = 0. Note that 
pi ,5 is hard to identify due to its small amplitude. (Inset) 
Prominent resonances of the spectrum F{uo) in dependence 
of N. The black solid lines are extracted from the Fourier 
transform of the fidelity F{t). The circles indicate the reso¬ 
nance positions obtained from our singlet theory (see text). 
Further, the excitation energies for the non-interacting case 
(g= 0 ) tuoi = ti — ei (dashed lines), hcoi = £3 — 61 (lowest 
dashed line), and fkou = 2 (e 2 — ei) (dashed dotted line) are 
shown. 
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the first trap state (see also the oblique magenta arrows 
in Fig. 0). Hence, it connects the inner part of the atomic 
cloud with the outer fraction establishing spatial coher¬ 
ence between them as we have seen in Fig. 0 Moreover, 
this mode induces, as we will see in the following, pop¬ 
ulation transfer between the inner and the outer atomic 
fraction. On the other hand, the harmonic oscillation 
can not be attributed to a single resonance, because no 
resonance shows up at its frequency. Nevertheless, we 
can understand its origin. We note that singlets (a^cq) 
with i = 5, 7,9,... to very high i are excited and there¬ 
fore present in F(u>). They correspond to oscillations 
between the numberstate \N, 0,0,...) and the states with 
| N — 1 , 0 ,..., 0 , 1 , 0 ,...), where a single particle is excited 
to the ith state and thus they oscillate, for g = 0 , with 
frequencies Wj = (e* — ei)//i with * = 5, 7, 9,.. (note that 
excitations into states with i even are symmetry forbid¬ 
den). These frequencies are shown in Figs. [5] as dashed 
lines. Since the difference between two neighboring reso¬ 
nances is approximatively constant ( 04+2 — uii) « 2 wo due 
to the equidistant spacing of the energy levels in the har¬ 
monic trap, all of these singlets are in phase again after 
a time T « tt/loo ~ 0.4 h/E*. Although this approxima¬ 
tion does not work perfectly at all times, because of the 
presence of the ionic potential, the harmonic oscillation is 
visible for the entire simulation time due to the contribu¬ 
tion of high energy modes which are only slightly affected 
by the presence of the ion. Nevertheless, the damping of 
the harmonic oscillation, visible in the trap energy (see 
also Fig. [3j green line), can be attributed to the induced 
slight dephasing. 

In the inset of Fig. 0 the resonance positions in depen¬ 
dence of the particle number N are shown. We observe 
that the peaks which can be associated to a coherence p\j 
shift to smaller energies for growing N. Our singlet the¬ 
ory is perfectly able to capture this behavior. Therefore, 
the shift of resonances can be understood as a mean- 
field phenomenon which effectively renormalizes the sin¬ 
glet resonances. Importantly, we numerically find that 
the property of equidistant spacing between the single¬ 
particle modes seems to be nearly untouched by the in¬ 
teraction and, as a consequence, the harmonic oscillation 
is essentially unaffected. 

The singlet theory explains most of the modes found in 
the Fourier spectrum F(ui). The low energy resonance at 
ui ss 15.5 E*/h, however, does not appear in the singlet 
spectrum. Moreover, a closer inspection of the inset of 
Fig. [5] shows that this mode has the tendency of an in¬ 
creasing energy as the number of bosons increases, which 
is the opposite behavior of the other resonances. We show 
in the following that this resonance can be attributed to 
correlations which are dynamically build-up during the 
temporal evolution. To do so we have first to understand 
the coupling between the singlets and the doublets. 

Towards that end, we proceed further with the anal¬ 
ysis of the dynamical evolution of the singlets. In Fig. 
[Gl the coherences P 13 , pn = )Cfc> 3 Pifc, and the occu¬ 
pations / 1 , / 2 , /t = J2k >3 fk are shown. As stated be¬ 


forehand, the coherence P 13 oscillates with the frequency 
of the ionic oscillation. From the time evolution of the 
ionic energy (see Fig. 0), one expects in addition the col¬ 
lapse and revival behavior. Here we see that for IV = 2 
(lower left panel), the absolute value of P 13 is nearly con¬ 
stant, while collapse and revival become clearly visible 
for IV = 10 (lower right panel). The coherence between 
the first bound and the trap states, Pit, shows the afore¬ 
mentioned rephasing behavior in the distinct peaks with 
T « 7 r/wo ~ 0.4 h/E* periodicity. The damping is also 
visible here for N = 10, but, importantly, the peaks are 
still very pronounced enabling the harmonic oscillation 
to be unaffected. Turning our attention to the dynamics 
of the occupations (Fig. 0 upper panels), we can iden¬ 
tify that the population of the bound states, thus the ion 
population f\ = f\ + / 2 , is about 80 — 90%, and there¬ 
fore only about 10 % — 20 % of the particles are emitted 
into the trap by the “ionization process”. Furthermore, 
a population transfer between the two bound states oc¬ 
curs (upper left panel). For larger N (upper right panel), 
also transfer to and from the trap population (green line) 
becomes visible with a rate equal to the frequency of the 
ionic oscillation. 

An inspection of Eq. (l20l) reveals that the population 
transfer between the bound state <t>\ and the trap states is 
mediated by coherences, primarily by P 13 , because it has 
the highest contribution in the Fourier spectrum F(ui) 
(see also Fig. 0, explaining its oscillation with the fre¬ 
quency of the ionic oscillation. In contrast, the dynamical 
evolution of fi has to be induced by the correlations T 22 , 
since all coherences in question vanish. Hence, the popu¬ 
lation transfer from the state <j>\ to the state <f >° would not 
take place in a mean-field scenario, thus it is a clear sig¬ 
nature of a genuine many-body effect. We also note that 
the frequency of this population transfer is w ps 8.2 E*/h 
for N = 2 which matches very well to the position of the 
additional low energy peak seen in the inset of Fig. 0 
Further, the larger N is, the faster the transfer process 
between the first bound and the trapped states happens 
which fits to the shift of the low energy resonance to 
larger energies visible in Fig. 0 We therefore show, in 
the next section, which doublets play an important role 
for this oscillation and to which process/mode it corre¬ 
sponds. Further, we will explain the origin of the collapse 
and revival of the ionic oscillation. 


C. Doublet Dynamics 

In order to understand the observed collapse and re¬ 
vival phenomena and the incoherent population transfer, 
we need to investigate the dynamical evolution of the 
doublets. For the sake of clarity, we classify the most 
important doublets (see definition in Sec. Ill B 21) into 
incoherent doublets gkkkk , gkqkq, and coherent doublets 
gkkqq, gkqkk • This separation is justified by the fact that 
the incoherent doublets are always real-valued quanti¬ 
ties, and therefore do not contribute to the dynamics 
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FIG. 6 . Time evolution of the singlets for N = 2 (left panels) and IV = 10 (right panels). In the upper panels, the occupations 
fk are shown: fa (black line), fa (red line), and fa (green line). Further, the absolute values of the most important coherences 
Pk q are plotted in the lower panel: jut as a black line, and P 13 as a gray shaded area (for better visibility). 


of the singlet occupations [see Eq. <H5D]- Moreover, 
the doublets have the properties gkqk'q' = dk’q'kq an d 

9kqk'q r — fjkqq'k' — Qqkk'q' — Qqkq'k'• Therefore, We 

can identify, noting that only correlations involving 
states up to k = 3 are contributing to the system 
dynamics, the most important classes of correlations: 
incoherent {31111,52222,53333,51212,51313} and coherent 
{31122,31133, 32233,31211,32122,31311,33133, 32322, 332331- 
In Fig. [3 we show for N = 2 (left panels) and N = 10 
(right panels) only those doublets of the incoherent 
(upper panels) and the coherent (lower panels) doublets 
which are considerably occupied. At first, we observe 
that, even though we start in an essentially uncorrelated 
state, very quickly a tremendous amount of correlations 
is created in the course of the temporal evolution. With 
the help of these correlations, we are now able to answer 
the remaining open questions concerning the dynamical 
evolution of the system. 

In order to understand the population transfer between 
the states (j>\ and (j> a , we have to identify the correla¬ 
tions g q kk'q' contributing to Fn and 1^2 ■ The coher¬ 
ent doublets 31122 ( 32211 ) is the only correlation in ques¬ 
tion which can be understood by inspecting Eq. (l22l) 
and Fig. [3 Comparing Figs. [3 and [3 (left panels), 
we see that they are build up and decay with the fre¬ 
quency of the population transfer between the state (f>\ 
to the state (fy- Since 31122 and 52211 which drive Fn and 
r 22 , respectively, are by 7r out of phase (31122 = 32211 ); 
they induce the population transfer between those two 
bound states. Thus, we can understand these coherent 


doublets as mediators for the population transfer. Fur¬ 
ther, we can now identify the low energy mode in the 
spectrum of Fig. [3 as a coherent oscillation between 
the numberstates \N, 0,...) and |IV — 2,2,0,...) which 
corresponds to a two-particle excitation. In the non¬ 
interacting case, this mode would therefore oscillate with 
a frequency wd = 2 (e 2 — ei )/H « 7AE*/H indicated by 
a dashed dotted line in the inset of Fig. [3 Here we can 
understand the following: First, this oscillation can not 
be seen in single-particle quantities like the one-body re¬ 
duced density matrix. Second, this is also the reason why 
we could not explain the low energy peak by our singlet 
theory. Third, any mean-field approach would not be 
able to capture the associated dynamics. For example, 
in a GP theory the population transfer from \N, 0,...) to 
|IV — 2,2, 0,...) would be not possible, because in a gen¬ 
eral GP state (X}fc c fc®!) iV l vac ) contributions from odd 
states are symmetry forbidden as long as only GP states 
with parity symmetry are considered. 

Finally, we would like to discuss possible reasons for 
the collapse and revival behavior of the ionic oscillation 
occurring for larger N. The impact of the doublets onto 
the singlet P 13 is contained in Even if many doublets 
contribute to Fi 3 , we can identify the responsible doublet 
by inspecting Fig. [3 (right panels). We observe that 
exactly at the times when the coherent doublets 31133 
are present, the ionic oscillation is strongly suppressed 
whereas when 31133 nearly vanishes the ionic oscillation 
reappears (compare Figs. [3 ® and[G|)- Therefore, the 
doublet 31133 is the main source of loss of the coherence 
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Pi3- Consequently, this doublet is responsible for the 
loss of spatial coherence between the inner and the outer 
density fraction during the dynamics, as it can be seen in 
Fig. 21 Nevertheless, the build-up of 31133 does not only 
act as a source of damping for the singlet coherences, 
since the coherences P 13 recur when the value of 31133 
is reduced again. Further, we note that 31133 (33311) is 
also present in Fn ^33) such that also the collapse and 
revival in the population transfer can be attributed to 
this coherent doublet. 


IV. CONVERGENCE ANALYSIS 

In the following, we briefly discuss the quality of our 
data and to which extend advanced tools as MCTDHB in 
the description of the dynamical behavior of such quan¬ 
tum systems are indeed necessary. To this aim, we in¬ 
spect the natural populations [see Eq. ©] which provide 
an assessment of the degree of fragmentation of the sys¬ 
tem [HI- Further, one can use them to systematically 
judge the convergence of our result to the exact solu¬ 
tion of the many-body quantum system [39,]. In Fig. 
0 the evolution of the natural populations exemplarily 
for TV = 10 particles is shown. We see that even if the 
initial state is nearly condensed, thus describable in the 
GP framework, depletion from the condensate state is 
created very quickly during the time evolution. There¬ 
fore, a GP description would be inaccurate and a multi- 
orbital description becomes unavoidable. Here we are 
using m = 6 single-particle orbitals. The contribution of 
the lowest natural orbital (smallest natural population) 
stays below 1% such that only natural orbitals with even 
smaller contribution are expected to be not taken into ac¬ 
count by truncating the Hilbert space. One could further 
speculate that in such cases a multi-orbital mean-field 
description might be sufficient. However, this is not the 
case and we stress that the non-steady behavior of the 
natural populations implies the necessity to go beyond 
a mean-field description [53]. Therefore, the utilization 
of MCTDHB is essential to capture the full correlated 
quantum dynamics of the system. 


V. CONCLUSIONS AND OUTLOOK 

We have investigated the dynamics of a trapped cloud 
of N interacting bosons after the sudden creation of a 
static ion. Thereby, most of the atoms are quickly cap¬ 
tured in the bound states of the atom-ion potential, while 
the remaining atomic fraction is emitted into the trap. 
We were able to understand the subsequent dynamics 
and the many-body spectrum in detail by means of an 
underlying singlet-doublet theory. 

Our atom-ion hybrid system exhibits, apart from the 
harmonic oscillation, that is, a density oscillation within 
the external harmonic confinement, an additional density 
oscillation which we were able to attribute to a coherent 


oscillation between a bound state and a trap state. This 
coherent oscillation results in the spatial coherence be¬ 
tween the inner and outer density fraction. It also allows 
for population transfer between bound and trap states 
such that the fraction of atoms in the bound states be¬ 
comes time-dependent. In contrast to the harmonic os¬ 
cillation, this ionic oscillation shows a particle number 
dependence. Furthermore, we showed that, even though 
the atoms are weakly interacting, strong correlations are 
build-up during the temporal evolution. These corre¬ 
lations lead to population transfer between the bound 
states and show up as a strong doublet resonance in the 
many-body spectrum. Besides this, the correlations in¬ 
duce a periodic suppression of coherence between the in¬ 
ner and the outer fraction which gives rise to collapse and 
revival of the ionic oscillation on longer time scales. 


The impact of the ion onto the bosonic cloud dynamics 
can be summarized as follows: Apart from the accumu¬ 
lation of atoms on both sides of the ion while depleting 
the cloud at the ion position, the ion induces a coher¬ 
ent oscillation between the bound states and the trap 
states. The oscillation essentially arises because of the 
additional scales present in the system. Indeed, this can 
be easily seen in the non-interacting case (3 = 0) and tak¬ 
ing the limit of small trapping frequency wo —> 0. In this 
case, the frequency of the ionic oscillation converges to 
the non-interacting single-particle energy e\ of the lower 
bound state (see Sec. Ill Al for definition) which, thus, sets 
effectively a lower bound to the frequency of the ionic os¬ 
cillation. Hence, although we have chosen a tight trap, 
which is numerically convenient to resolve the dynamics 
in a finite system, the ionic oscillation would be present 
with a frequency on the same order of magnitude for a 
weak confinement, too. 


The present work can be viewed as a further step to¬ 
wards the simulation of the dynamics of the hybrid atom- 
ion system in the ultracold regime. The fact that quasi 
one-dimensional quantum Bose gases are routinely re¬ 
alized in several laboratories and given the recent ad¬ 
vances in optical trapping of ions, puts the experimental 
realization of the hybrid system considered here within 
reach. Furthermore, since the multilayer extension of 
our method (ML-MCTDHB) is especially designed for 
the simulation of mixtures, we plan to investigate the 
impact of the ionic motion onto the ground state and the 
dynamical evolution of the bosonic ensemble in the nearer 
future. This will be done by treating the ion quantum 
mechanically as well. Due to the attractive inter-particle 
interaction, these future studies can reveal the dynami¬ 
cal formation of molecular ions as predicted in Ref. [23| . 
Furthermore, it would be of interest to study the energy 
transfer between ionic and atomic degrees of freedom in 
order to understand, for example, sympathetic cooling of 
the ion in the atomic cloud in more detail. 
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FIG. 7. Time evolution of the most important doublets for N = 2 (left panels) and N = 10 (right panels). In the upper panels, 
the incoherent doublets <71111 (black solid line), <72222 (red dashed line), and <73333 (blue line with crosses) are shown while the 
absolute values of coherent doublets <71221 (black solid line), <71122 (red dashed line), and <71133 (blue line with crosses) are shown 
in the two lower panels. 
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Appendix A: Derivation of Clusters from Density 
Matrices 

The M-particle clusters can be calculated via the M- 
particle reduced density matrices. Let us start with 
the time-dependent one-particle reduced density matrix 
which describes the spatial density and coherence of the 
system. We express it in terms of the natural orbitals 
and the natural populations A j as 

p(x,x',t) = ^ \ j (t)$*(x,t)$ j (x',t). (Al) 

j 

If we now expand the natural orbitals in an arbitrary 
single-particle basis <j)j{x) 

=Y J C j k(t)M*) ( A2 ) 

k 


with 

Cjk{t)= j <t>* k {x)$j{x,t)dx (A3) 

and express also p(x,x',t) in this basis 

p(x, x\ i) = £ (f>* k {x)(j) q {x')(ala q ), (A4) 

k,q 

we can identify the single-particle cluster as 

( A5 ) 

3 

Here aj, and a k are the creation an annihilation operators 
of the single particle basis <f>k(x). In the same manner, we 
can proceed in order to derive the two-particle clusters. 
By defining the expansion coefficients as 

d kq( t ) = J <t>* k (xi)(j)* q (x2)®j{x 1 ,X2,t)dx 1 dx2 7 (A6) 

where now (xi , X2 , t) are the geminals, that is, the 
eigenfunctions of the two-body density matrix, we obtain 
the following expression for the two-particle clusters 

{a\a\a q ,a k ,) = ( A7 ) 

3 
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We obtain 


ih^5(a\aj) = ^\fl q 5 ki - e°^S qj 

kq 


H" 2 ^ ^ ^kjk'g{0'jQ'k/) 0 
k' 

2 ^ ^ Vqk'ikiO'fcrQ'j )0 
k' 

“ 1 “ {y^kkki “ 1 “ ^ ^ Vq'kki{Q' a 'Q J k)o')3jq 
q' 

— (Viqqq + El ^j'g<?g' (QoQ Q ')o)^fc 
q' 

+ Vkgqiialaj )0 - V^fcfc^ala^o <5(aJ.a g ) 

+ f« (B2) 


where e°- = e^j + 2 Efc, V ki j q {a\a q )o and T R0 coincides 
with the definitions of F R . but with the mean values of 
the singlets inserted, whereas F^ is defined as 


FIG. 8 . Natural populations A, (£) exemplarily for IV = 10 
particles on a logarithmic scale. Note that the A j are by 
definition ordered such that Ai > A 2 > ... > Am for all times. 


Appendix B: Linearized Singlet Dynamics 


F — F 

L *j — ij 


■^BO 


q 


(B3) 


Here we used Eq. (ED for the correlations T R . This is an 
important step because in case the doublets are negligible 
it decouples the singlets dynamics from any other time- 
dependent driving such that the many-body system can 
be described by singlets only. Finally, Eq. (IB2I) can be 
rewritten in matrix form as 

ifi-S(t) = MS(t) + f It). (B4) 

df 


In order to derive the energy spectrum of the singlets, 
we need to linearize the equations of motion (ED and 
(ED- Even though such an approximation can not be 
able to describe the full dynamics, we should be able to 
predict the spectrum of the singlets rather accurately. To 
this end, we write the singlets as 


such that the homogeneous solution can be expanded in 
terms of the eigenvalues v n and the left(L) or right (R) 
eigenvectors S RR of the matrix M 

MS r = i/„S R , (B5) 

and the imhomogeneous part can be found by means of 
the variation of the constant method leading to the full 
solution 


(«K->(*) = (ajaj) 0 + S(alaj){t), 


(Bl) 


S(*)=E 

n 


C-n 


ft 1 

/ — e lUnt /ri S(ff (t')dt' 
Jo ^ 


qR —iVnt/h 
&n e 


(B 6 ) 

The coefficients c n determine how strong the mode S R 
is excited and can be obtained by projecting the above 
solution onto the eigenvector basis, resulting in 


where (a|a -)0 is the time-independent mean value of 
(ajci - )(i) and 8{a\a^){t) is a small time-dependent fluctu¬ 
ation. By neglecting the terms quadratic in 5{a\a^), we 
get a system of linear differential equations for 8{a\a ■). 


ft I ; 

^ - / — e iVnt / ri S(ff (t')dt'. (B7) 

Jo ™ 

In summary, we see that the linearized singlet equa¬ 
tions provide us with the singlet resonances v n and the 
singlet inodes S R . We can derive these from the MCT- 
DHB solution in the following way: Via Eq. (IA5I) . we can 
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extract the dynamics of the singlets which we can use to 
derive the matrix M, and therefore the solution S(f). By 
diagonalizing M, we obtain the singlet energies and the 
dominant singlet in the eigenmodes (largest entry in the 
vector) shown in Fig. [5] Further, we can extract the 


oscillator strength |c„| 2 by Eq. (lB7l) which we can use 
as a measure for the relative height of the peaks in the 
many-body Fourier spectrum of F(t) just by scaling them 
to the global maximum of the spectrum. Note that, Eq. 
(IB7I) should be evaluated at small t in order to allow the 
linear approximation. 
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